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Abstract This paper extends the energy-based version of the stochastic hnearization method, 
known for classical nonlinear systems, to open quantum systems with canonically commuting dy- 
namic variables governed by quantum stochastic differential equations with non-quadratic Hamil- 
tonians. The linearization proceeds by approximating the actual Hamiltonian of the quantum 
system by a quadratic function of its observables which corresponds to the Hamiltonian of a quan- 
tum harmonic oscillator. This approximation is carried out in a mean square optimal sense with 
respect to a Gaussian reference quantum state and leads to a self-consistent linearization proce- 
dure where the mean vector and quantum covariance matrix of the system observables evolve in 
time according to the effective linear dynamics. We demonstrate the proposed Hamiltonian-based 
Gaussian linearization for the quantum Duffing oscillator whose Hamiltonian is a quadro-quartic 
polynomial of the momentum and position operators. The results of the paper are applicable to 
the design of suboptimal controllers and filters for nonlinear quantum systems. 

1 Introduction 

A wide class of models for open quantum systems [[31 [8l, that is, quantum-mechanical ob- 
jects interacting with the environment, is provided by dynamical systems whose state vari- 
ables are canonically commuting self-adjoint operators on a Hilbert space. In the Heisen- 
berg picture, these system observables evolve in time according to quantum stochastic 
differential equations (QSDEs) ll23ll . Such QSDEs, which are dual to quantum master 
equations for density operators in the Schrodinger picture flSl Chapter 3], are driven by 
a quantum Wiener process to take into account the coupling between the environment 
(regarded as a memoryless heat bath of quantum harmonic oscillators) and the internal 
dynamics which the system would have in isolation from the surroundings. These inter- 
nal dynamics are completely specified by the system Hamiltonian, which is a self-adjoint 
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operator on the underlying Hilbert space, usually representable as a function of the system 
observables. 

In particular, quadratic system Hamiltonians correspond to quantum harmonic oscil- 
lators whose behaviour lends itself to complete analysis due to linearity of the resulting 
QSDEs in contrast to the general nonlinear case. Linear open quantum systems are being 
actively researched to develop quantum analogues of classical control schemes, including 
the Tioo, risk-sensitive and linear quadratic Gaussian control approaches (see, for exam- 
ple, [|71 Il2l [121 [I6l |22l [30l [m and references therein). Such models are also employed in 
quantum optics which is considered to be one of possible platforms for implementing the 
quantum computer [I2T1 Section 7.4]. 

The present paper is aimed at a quantum-mechanical version of the stochastic lin- 
earization (SL) technique whose origins date back to [|2l [51 [141 (see also fTl 161, l27l and 
references therein). SL is concerned with a classical SDE whose drift term is a nonlin- 
ear function of the state vector. The principal idea of SL is to approximate the drift by an 
affine function of the state variables whose coefficients are computed using a mean square 
criterion with respect to a probability distribution. This reference distribution, which is 
intended to mimic the actual probability distribution of the state vector, is usually chosen 
to be Gaussian, although non-Gaussian approximations (such as, for example, in [j6ll27ll) 
are also utilized. The Gaussian reference measure leads to an effective linear SDE which 
approximates the actual nonlinear dynamics. A salient feature of this SDE is that its co- 
efficients depend nonlinearly (through integral operators with Gaussian kernels) on the 
mean value and covariance matrix of the state vector, which are in turn governed by linear 
ordinary differential equations (ODEs) (including the Lyapunov ODE for the covariance 
matrix) involving those coefficients!^ This provides a self-consistent procedure for lin- 
earizing the dynamics. 

An alternative energy-based version (3T\ of the SL technique was aimed originally at 
structural engineering problems of random vibrations with a potential nonlinear restoring 
force. Rather than directly linearizing the nonlinearity, this approach employs a mean 
square criterion in order to approximate the force potential by a quadratic function of the 
displacement vector (corresponding to an ideal spring). It is this variant of the classical 
SL that is particularly suitable for our purposes. We adapt it to the quantum-mechanical 
setting by solving the problem of minimizing the mean square deviation between the ac- 
tual non-quadratic system Hamiltonian of the quantum system and a general quadratic 
function of its observables. The solution involves the second and higher-order mixed 
moments, which, in a Gaussian quantum state [|24| (see also [[8l pp. 1 18-122]), are com- 
pletely specified by the mean vector and the quantum covariance matrix of the system 
observables through Wick's theorem [17, p. 122]. 

For a class of open quantum systems, whose coupling with the external heat bath 
variables in the total Hamiltonian is bilinear, the quadratic approximation of the system 
Hamiltonian leads to a linear QSDE of an open quantum harmonic oscillator which is 

'This resembles the McKean-Vlasov SDE (from the kinetic theory of plasma) whose drift depends on 
the probability density function of the state vector propagated by the Fokker-Planck-Kolmogorov equation 
associated with the SDE, thus leading to a nonlinear parabolic partial differential equation [19). 
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amenable to comprehensive analysis. In particular, similarly to the classical linear systems 
[[T5ll . the quantum covariance matrix of the observables of this effective oscillator satisfies 
a Lyapunov ODE. Moreover, such linearization respects the physical realizability (PR) 
conditions [|T3l Theorem 3.4 on p. 1790], which makes it suitable for coherent quantum 
control II22II . We demonstrate the approach for the quantum Duffing oscillator flU |25l 
with a quadro-quartic Hamiltonian. The proposed Hamiltonian-based quantum Gaussian 
linearization technique is applicable to the development of suboptimal controllers and 
filters for nonlinear quantum systems since it offers a recipe to deal with the "curse of 
dimensionality" of the information state, similar to projective quantum filtering [|29l . 

The paper is organised as follows. Section [2] specifies the class of quadratic Hamilto- 
nians. Section |3] describes the approximation of an arbitrary Hamiltonian by a quadratic 
Hamiltonian, optimal in the mean square sense. Section |4] specializes the computations 
for a Gaussian quantum state. Section [5] describes the Hamiltonian-based self-consistent 
linearization for open quantum systems. Section |6] demonstrates the hnearization proce- 
dure for the quantum Duffing oscillator. Section |7] provides concluding remarks. Long 
proofs and subsidiary material are given in Appendices. 

2 Quadratic Hamiltonians in canonically commuting vari- 
ables 

Suppose xi, . . . ,Xn are quantum observables (that is, self-adjoint operators on an under- 
lying separable Hilbert space "H with an inner product {ip \ ^)^which satisfy canonical 
commutation relations (CCRs) 

[xj, Xk] = i9jkl, 1 ^ j,k ^n. (1) 

Here, i := is the imaginary unit, [A, B] := AB—BA is the commutator of operators, 
and 9 := {6jk)i^j^k^n is a real antisymmetric CCR matrix of order n (the space of such 
matrices is denoted by A„). Also, X denotes the identity operator which carries out the 
ampliation of entries of the matrix 6 to the space of linear operators on 1-L and will be 
omitted for brevity, so that ([T]) can be written in a vector-matrix form as 

[x, x^] := ([xj, Xk])^^j^k^^ = iO, (2) 

where the observables are assembled into a vector x := {xj)i<^j<^n- Unless indicated other- 
wise, vectors are organised as columns. The transpose (■)^ applies to vectors and matrices 
with operator- valued entries as if the latter were scalars. In particular, the CCRs hold for 
self-adjoint operators which are representable as linear combinations of annihilation and 
creation operators ai, . . . ,au and a\, . . . ,al, where u := n/2 and n is assumed to be even. 
Such are, for example, the quantum-mechanical position and momentum operators q and 

^To avoid confusion, we use a different notation (•, •) for other inner products, for example, the Frobe- 
nius inner product of matrices. 
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p := 



—idg with [q,p] = i and CCR matrix 



which spans the space A2. The associated annihilation and creation operators [|28l pp. 
90-91] 

a := {q + ip)/\/2, := {q — ip) / V2 (4) 

satisfy [a, a^] = 1. Now, with a scalar a G M, a vector b := ibj)i<^j<^n ^ IR" and a real 
symmetric matrix R := (rjfc)i<jj of order n (the space of such matrices is denoted by 
§„), we associate a self-adjoint operator 

Ha,b,R ■.= a + b^x + x^Rx/2 = a + ^(bk + - ^ Tj^x^Xk (5) 

i=i i=i 

on the Hilbert space 1-L. The operator Ha^h R, which is parameterized linearly by the 
triple (a, b, R) E M x M" x §„, is the Hamiltonian of a quantum harmonic oscillator 
with state variables xi, . . . , x„. Although the constant term a in ([5]) has no influence on 
the system dynamics, it is retained to preserve the generality of Ha,b,R as a quadratic 
polynomial of the system observables with real coefficients. If the system is isolated from 
the environment, its Heisenberg dynamics are described by the ODEs 



xt = i[Ha,b,R, xe] = ^( bj[xj, xe] + ^^ rjk[xjXk, Xi]^ 

j=l j,k=l 
n ^ n n n 

= -Yl ^^^^^ ~ 2 S rjkiOkeXj + djiXk) = Oij (^bj + ^ rjkXk) , (6) 

j=l j,k=l j=l k=l 

where the commutator identity [AB, C] = A[B, C] + [A, C]B from ['20'. Eq. (3.50) on p. 
38] is applied to the triple Xj, Xk, xi and use is made of the CCRs from ([!]) along with the 
antisymmetry of 0. In a vector-matrix form, (|6]) can be written as 

X = i[Ha,b,R, x] = e{b + Rx). (7) 

If i? :^ 0, the spectrum of the matrix QR is purely imaginary, and the system is neutrally 
stable. In general, R is not necessarily positive definite. For example, quantum amplifiers 
flSJ, used as active elements in quantum optics, are modelled as inverted oscillators with 
i? -< 0. If i? is nonsingular, the effect of b reduces to a constant shift R^^b in x, so that ^ 
can be written in terms oiy := x + R~^b as y = QRy. If the system has a non-quadratic 
Hamiltonian H, then, in contrast to (|7]), the right-hand side i[H, x] := {i[H, Xj])i^j^n of 
the Heisenberg dynamics is not affine in the system observables. In this case, i[H, x] can, 
in principle, be approximated by 0(6 + Rx) so as to minimize a mean square deviation 
between the vectors: 

E{Al^FAb,R) min, b e M", R G §„. (8) 
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Here, EA := Tr(pA) denotes the quantum expectation of a linear operator A on the 
underlying Hilbert space "H with respect to a density operator p (a positive semi-definite 
self-adjoint operator on V. with unit trace Trp = 1) which specifies the quantum state [|23l 
p. 51]. Also, 

Ab,R:=i[H,x]-e{b + Rx) (9) 

is a vector of "residuals" which depends affinely on b and R, and F is a complex positive 
definite Hermitian matrix of order n. Such linearization of i[H, x] can be regarded as a 
quantum version of the weighted least squares from the classical linear regression analysis 
j26l. The weight matrix F, which governs the minimization problem ([8])-(|9]), influences 
the optimal values of b and R, and its particular choice requires additional consideration. 
We will therefore take a different approach to hnearizing the system dynamics, through a 
quadratic approximation of the Hamiltonian itself, as a quantum counterpart to the energy- 
based variant of SL given in [|32|| . 



3 Mean square optimal quadratic approximation of Hamil- 
tonians 

Consider the mean square optimal approximation of the system Hamiltonian H hy a 
quadratic Hamiltonian Ha^b^R from ([5]): 

Q{a, f3, R) := E{{H - H,^b,R?) = - K,p,R?) min . (10) 

Here, we have introduced a different parameterization of the quadratic Hamiltonian 

Ha^hfi = a + l3^i + =: K,^^n (1 D 

through "centering" the observables of the system in the reference quantum state: 

^■.= iQi^j^n:=x-Ex, r,:=H-EH. (12) 

The new parameters a := a + b'^Ex + (Ex)^REx/2 — EH and b + REx are bijectively 
related to the old ones a and b from (|5]), with the matrix R remaining the same. To 
compute the optimal values of a, (3, R which minimize Q{a, /3, R), we will use the real 
parts of the following mixed central moments of the actual Hamiltonian H and the system 
observables xi, . . . , Xn- 

ej := ReE{r]Q, := ReE{r]^j^k), (^jk ■= ReE{^j^k), (13) 

TjM := ReE{^j^k^(,), (pjkim ■= ReE{^j^k^e^rn)- (14) 
Several remarks are in order on the matrices F := {'jjk)ifS,j,k^n and S := {crjk)i^j,ki^ 



and the tensors T := (r,H)is£j,fc,£ssn and $ := {(pjkem)i<^j,k,e,m^n defined by (13)-((T4|). 



Since the system observables satisfy CCRs, then in view of Lemma [2] from Appendix |A[ 



5 



the matrix T = ReE(r7^^^) is symmetric and the tensor T is totally symmetric. In what 
follows, T is identified with a linear operator acting from §„ to M" as 

n 

T{R) := ( y] Tjkerke) , R := {rke)i^k,eiin e §n- (15) 

The matrix S = ReE(^^^) is symmetric and positive semi-definite as the real part of 
the quantum covariance matrix of observables. However, the tensor $ from ( 14) is only 
guaranteed to be symmetric with respect to reversing the order of its subscripts. Indeed, 
similarly to ( A3 ) from Appendix [Af 



with (■) the complex conjugate, and hence, the quantum expectations on the opposite 
sides have equal real parts, that is, (fjkem = ^mikj- The fact that (^jkim is, in general, 
not invariant even under transpositions of its neighbouring subscripts follows from the 
identities 

'^jMm — ^kjtm = 'R'GE{[^j,^k]^£^m) = ^G{i9jkE{^e^rn)) 

= Re{i9jk{aem + idtm/'^)) = —OjkOhn/'^ = ^jklm — '^jkml (16) 

and a similar relationship 

^jkem — y^jekm = —Ojmdkl/'^-, (17) 

which are established by using the CCRs ([T]) and the definitions ( 13 ), ( 14). Now, consider 
a partial symmetrization ^ := {ip jMm)i^j,k/,m^n of $ whose entries are defined by 

i^jkim '■= i'^jkim + ^jkmi + '^kjim + fkjme 

+'^emjk + '^mljk + ^Imkj + ^mlkj)!'^, (18) 

which involves only eight of the 24 possible permutations of the subscripts j, /c, m. We 
will identify ^' with a self-adjoint operator on the Hilbert space §„ (with the Frobenius 
inner product of matrices (X, Y) := Tt(XY) inherited from M"^") defined by 

n 

^(i?) := ( ipjkemremj , R ■= {rem)i^e,nis:n G Sn- (19) 



'm=l 



The significance of \E' as the partial symmetrization of $ (with the latter being regarded 
as a linear operator on M"^", defined similarly) is that 

EHeRO') = {R, HR)) = {R, ^m, R e Sn, (20) 

which, in fact, can be used as an equivalent definition of \1'. Since the quantum expectation 
of a squared observable is always nonnegative, pO]) implies that the operator ^ is positive 
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semi-definite (\& ^ 0). By using the identity E{AB) = E{BA) for observables A, B, it 
follows that the mean square criterion < [T0] ) is a convex quadratic function: 

Q{a, P, R) = E(r/2) - 2ReE{r]h^,p,R) + E{h'). (21) 



In view of ([TT[)-([T3]), the second term on the right-hand side of ( pT[ ) does not depend on 
a and is linear with respect to (3 and R: 



ReE(r/V/3,R) = e^/3 + (r,i?)/2. 



(22) 



By a similar reasoning, ( 1 1 )-([T4]) imply that the rightmost term in (21 1 is a positive semi- 
definite quadratic form 



(C,n(C)>, C:^ 



a 
R 



which is specified by a self-adjoint operator 11 on the Hilbert space M x 
the inherited inner product (C, (') ■= + Z?^/?' + {R, R')) as 



n(C) := 



a + (S,i?)/2 
^l3 + T{R)/2 
Sa/2 + 7^(/3)/2 + ^'(i?)/4 



1 (5],-)/2 

s r/2 
s/2 rV2 



(23) 
X S„ (with 

C- (24) 



Here, use is made of the operators T and \l/ from ( 15 1, ( [l9| ), and the adjoint operator : 
M" §„ maps a vector (3 := (/3j)i^jXn to a matrix T^{(3) := ( Zlfc £=1 ^jh/^jOk^ 



The symbolic matrix representation of 11 in ( 24 ) can be identified with the real part of 
the matrix of second moments of the triple (1, ,^,^^/2). By a generalized version of 
the Schur complement condition of positive definiteness [l9l Theorems 7.7.6, 7.7.7 on pp. 
472-474], the invertibility of the positive semi-definite operator 11 is equivalent to S 
and G y 0, where G is a self-adjoint operator on §„ defined by 



' 1 ' 


-1 




s 




T 



^-S(S,-) -T'^S^^r. (25) 



Up to a factor of 4, the operator G is the Schur complement of the block 



1 
s 



in (24). 



Theorem 1 Suppose the operator IT, defined by (24), is invertible. Then the optimal 



values of the parameters a, (3, R, which minimize the fixnction Q{a,(3,R) in (21), are 



computed in terms of (13), (14) and (25) as 



R(, 



-(S,i?o)/2, 
S-i(e-r(i?o)/2), 
2G-\T -T\T.-^t)). 



(26) 
(27) 
(28) 
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Proof. The optimal values of a, /3, R are obtained by equating the Frechet derivatives of 



the function Q{a, (3, R) in (21 1 to zero. In view of (22) and (23 1, this leads to the system 
of linear equations 



d^Q = 2a+(S,i?) = 0, 
dpQ = 2S/3 + r(i?) - 2e = 0, 
OrQ = *(/2)/2 + r^(/3) + «S 



0, 



(29) 
(30) 
(31) 



which correspond to the normal equations of the least squares method in the linear regres- 
sion analysis [|26ll . If the operator 11 ^ in ([24]) is invertible (and hence, 11 >- 0), then 



the quadratic function Q in (21 1 is strictly convex and the system of equations (29)-(31 ) 
has a unique solution. Now, (26) and ([27]) follow directly from (29) and (30), while (28) 
is established by their substitution into (ST) and using the invertibility of the operator G 



from (25 ) which is secured by the condition 11 ^ 0. 



Since the vector e, the matrices T and S and the tensors T and ^, which are associ- 
ated with the mixed central moments of the Hamiltonian H and the system observables 
Xi, . . . , Xn, depend on the quantum state, then so also do the parameters a^, R^ of the 
optimal quadratic approximation of the Hamiltonian. 



4 Approximating the Hamiltonian in a Gaussian quan- 
tum state 

The system is said to be in a Gaussian quantum state [|24ll . if the quantum covariance 



function of the centered vector ^ of system observables from ( 12) is given by 
Here, 

S := {sjk)i^j,k^n := m^^) = S + ze/2, (33) 

is the quantum covariance matrix, which is a complex positive semi-definite Hermitian 
matrix, with 6 and S defined by ([2]), (13), and use is made of the antisymmetry of 6. 



By applying Wick's theorem pTl p. 122], which is a quantum counterpart to Isserlis' 
theorem [fTOl on the mixed central moments of evenly many jointly Gaussian classical 
random variables in terms of their covariances (see also [ITTl Theorem 1.28 on pp. 11- 
12]), it follows that, in the Gaussian quantum state, 

r 

mn X • • • X O^J = E n ^^-^..-.i^.. ■ (34) 



Here, the sum of products of the quantum covariances from (33 ) extends over a class Vr of 



(2r— 1)!! permutations {ki, . . . , k2r) of the integers 1, . . . , 2r which satisfy k2e-i < k2e for 
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every 1 ^ i ^ r and ki < < . . . < k2r-3 < A;2r-i- Such permutations will be referred 
to as regular. There is a one-to-one correspondence between the regular permutations and 
all possible partitions {{ki, ^2}, . . . , {A;2r-i, fc2r}} of the set {1, . . . , 2r} into two-element 
subsets. Thus, ([34]) allows any mixed moment of even order 2r to be computed for the 
observables ^1, . . . in a Gaussian quantum state in terms of the matrix S, whereas 
all the moments of odd orders in such a state are zero; see Appendix |B] In particular, 



the tensor T of third order moments from ( [T4| ) vanishes, while application of ( [34] ) to the 
fourth order mixed moments yields 

E(^j'Cfc^^^m) = SjkSlm + SjlSkm + SjmSkl, (35) 

cf. a similar relation for the annihilation and creation operators (]4]) in [[8l Eq. (4.4.121) on 



p. 122]. Hence, the real parts of the fourth order moments in ( 14) take the form 



which is in accordance with the more general relationships (16), ( 17). However, it will 



be more convenient to compute the partial symmetrization \E' of $ by applying (35 ) to the 



equivalent definition of ^ in (20) rather than using the entry wise representation (36) 



Lemma 1 Suppose the system is in a Gaussian quantum state. Then the operator ^, 



defined by {18)-{20) in terms of the tensor from ([i4[), takes the form 



where K is a positive semi-definite self-adjoint operator on the space §„, defined by 

K{R) := + ei?e/4, 



(37) 



(38) 



with and S defined by ([2]) and (13). If the quantum covariance matrix S from (33) is 
nonsingular, then K y 0. 

We prove Lemma[T]in AppendixJC] In a particular case, when the system observables 
commute with each other, that is, 6 = 0, Lemma [T] reduces to the well-known result on 
the second moment of a quadratic form in jointly Gaussian classical random variables 
[[T8l Lemma 2.3 on p. 204]. If 6 7^ 0, the noncommutative quantum nature of the 



system observables enters (38) through the additional term 0i?0/4 which makes K a 
special self-adjoint operator of grade two [30, Section 7]. The above discussion allows 
Theorem [T] to be concretized for the Gaussian quantum case as follows. 



Theorem 2 Suppose the mean square deviation Q(a, j3, R) in {21) is associated with a 



Gaussian quantum state, and the quantum covariance matrix S in {33) is nonsingular. 
Then the optimal values of a, (5, R in (]26])-([2^ are computed in terms of the mixed 



central moments {13 ) and the associated positive definite operator K from {38) as 



-(S,i?,)/2, (3, 



Ro 



(39) 
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Proof. The expressions ( [39| ) are obtained from (|26|)-(|28)) by noting that, in the Gaussian 
quantum state, T = and, in view of Lemma [T] the operator (25 1 takes the form G = 
^ — •) = 2K, where the invertibility of K is ensured by the assumption that S' 0. 



Although the operator in (38) is completely specified by the matrices S and 6, 



the optimal parameters /3o and in (39) will also depend on the mean value Ea; in the 
Gaussian reference state through the vector e and the matrix T. The computation of the 



inverse operator K~^, which is required for (39), is described in Appendix |D] where it is 
also shown that the condition S y ensures the positiveness of K^^ with respect to the 
convex cone of real positive semi-definite symmetric matrices of order n in the sense 
that 



(40) 



5 Self-consistent quantum Gaussian linearization 

Suppose the quantum system interacts with the external heat bath so that the n-dimensional 
vector Xt of its observables at time t is governed by a QSDE 

dXt = {i[H, Xt] - BJB^e~^Xt/2)dt + BdWt, (41) 

where H is the system Hamiltonian discussed previously, and the CCR matrix 6 of the 
system observables is assumed to be nonsingular. This corresponds to a bilinear coupling 
between the open quantum system and the bath variables in the total Hamiltonian as 
quantified by a constant matrix B E M"^™; see, for example, [|7l[T3l for details. Also, 
Wt is an m-dimensional quantum Wiener process (with m even) which represents the 
influence of the environment on the system. The entries of Wt are self-adjoint operators 
on a boson Fock space [|23l with the quantum Ito table 



dWtdW;" = Qdt, n:=Ira + iJ/2, J = J®/„/2, (42) 

where Ir denotes the identity matrix of order r, the matrix J is given by (|3]), and ® is 
the Kronecker product of matrices, so that J is the CCR matrix of Wt in the sense that 
[dWti dW'^] = iJdt. Consider the first two moments of the system observables: 

fit:=EXt, St:=Re5i, St := E{^t^J), ^t ■= Xt - fif (43) 

Although the quantum state of the system is not necessarily Gaussian, fxt and can 
be used to compute the parameters of the mean square optimal quadratic approximation 
Hao,i3^,Ro of the actual Hamiltonian H through Theorem [2] as if the state were Gaussian. 



In this Gaussian reference quantum state, the term i[H, Xt] in (41 ) can be approximated 
as 

i[H, Xt] ^ i[H^Mo, Xt] = i[(3j^t + 6] = e(/3o + RM)- (44) 
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Formal substitution of ( [44| ) into the right-hand side of ( [41] ) yields a linear approximation 
for this QSDE which splits into an approximate ODE for the mean jit and an approximate 



QSDE for the centered vector of system observables from (43 ): 

[it ^ Q(3, - BJB^&-'iit/2, d6 ^ A6dt + BdWt. 
Here, the external field is assumed to be in the vacuum state [|23l , and the matrix 

At := ei?o - BJB^e-^/2 



(45) 



(46) 



depends deterministically on time t through which is completely specified by jit, ac- 
cording to Theorem[2| Upon averaging, the quantum Ito differential d{C,tC,t') = (d.^,)^,^ + 
^td^'t + (d'Ct)d'Ct^, combined with (45) and (42) according to [1231 Proposition 25.26 on 
pp. 202-203], leads to an approximate Lyapunov ODE St ^ AtSt + StAj + BQB'^ for 



the quantum covariance matrix St from ([43]). If this approximation is regarded as an exact 
equation 



St = AtSt + StA'+BnB' 



(47) 



its solution satisfies St = Tjt + i&/2 y for all t ^ 0, provided So := Sq + iQ/2 >- 0. 
Here, we have used the positive semi-definiteness of the matrix f2 from ( |42| ) and the 
monotinicity of the transition operator of the Lyapunov ODE, with the fact that At from 



(46) satisfies 



AtQ + QAj + BJB^ 



0. 



(48) 



The latter property, which is equivalent to the preservation of the CCR matrix in time, is 
one of the physical realizability (PR) conditions f\3[ Theorem 3.4 on p. 1790] describing 
the dynamic equivalence of the system to an open quantum harmonic oscillator. Now, if 
the first of the equations ( [45] ) is also considered as an exact ODE, then its combination 
with the real part of ([47]) yields 



e/3o - BJB^Q-^jit/2, tt = AtEt + EtAj + BB^ 



(49) 



In conjunction with (46), the ODEs ( [49] ) provide a self-consistent set of nonlinear equa- 
tions for finding jit and as functions of time, with the nonlinearity coming from the 



dependence of /3<^, in ([39]) on jit, S,. Therefore, although the above ODEs result 
from an ad hoc approximation, this quantum Gaussian linearization procedure generates 
a faithful quantum covariance matrix for the system observables and respects the PR con- 



ditions in the sense of (48). A time invariant version of the procedure is obtained by 



equating the right-hand sides of (49) to zero and considering admissible solutions ji, S of 
the corresponding algebraic equations 



Ql3^- BJB'e-'ji/2 = 0, AS + SA 



T 



BB^ 



0, 



for which the matrix 



A := ei?o - BJB^Q-^/2, 



(50) 



(51) 
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obtained from ( |46l ), is Hurwitz. If the matrix B is of full row rank, then, in view of i7 :^ 0, 
the corresponding solution S of the algebraic Lyapunov equation AS + SA^ + BVlB^ = 
as a steady- state version of ( [47| ) (assuming A Hurwitz) is nonsingular, thus ensuring the 
invertibility of the operator K, which is essential for the quantum Gaussian linearization 
through Lemma [T] and Theorem |2] Since the functions /3o and can, in general, be of 
complicated nature, their presence makes (50)-([5T|) a system of nonlinear vector-matrix 
algebraic equations for which the existence/uniqueness of admissible solutions yU, S is a 
nontrivial problem. Nevertheless, the example in the next section demonstrates tractabil- 
ity of this problem for a class of anharmonic oscillators. 



6 Application to the quantum Duffing oscillator 

Consider the quantum Duffing oscillator flU |25l of unit mass on the real line with the 
Hamiltonian 



X 



q 
p 



(52) 



Here, the vector x of system observables is formed by the position and momentum opera- 
tors q and p from Section [2] with the CCR matrix := J given by Q. The quadro-quartic 
polynomial V{q) describes the potential energy, where luq is the harmonic frequency and 
the coefficient / "weights" the quartic part, which is responsible for the anharmonicity of 
the oscillator. Application of the relation = V'{q)[q,p] = iV'{q), which follows 

from the commutator identity of [I20l Eq. (3.51) on p. 39]) and the CCR [q,p] = i, leads 
to 



i [H, x] 



P 



containing a cubic term. We will now demonstrate the Hamiltonian-based quantum Gaus- 
sian linearization technique of Sections |4] and [5] for the quantum Duffing oscillator. Ac- 
cording to ( E16[ ) and ( E17[ ) of Appendix |E} whose derivation is based on Theorem [2] and 
repeated use of Wick's theorem, the parameters of the mean square optimal quadratic ap- 



proximation of the Hamiltonian (52 1 in a Gaussian reference quantum state are computed 
as 



l3o = Nfi, N :-- 



a;2 + 4/(^2 + Sail) 
1 



1 



(53) 



Here, /i := Ea; 



Ep 



, with K := Eg and cth := E((g — k)^) the mean and variance 



of the position operator. Suppose the open quantum Duffing oscillator is governed by the 
QSDE (41 1, with B E M^^™ of full row rank. Since the matrix J spans the space A2, 
such that 



there exists a real 



BJB^ 



0J. 



(54) 
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Hence, BJB^Q ^ = 0/2, and the time invariant version of the quantum Gaussian lin 
earization procedure ( [49] ) reduces to the algebraic equations 



where 



(JA^ - 0/2/2)/i = 0, AS 



C '.— (CjA:)l^j,fcs£2 



0, 



BB^ 



satisfies C + icpJ /2 = BVtB >- 0, and the matrix A from (46) takes the form 



A := 3R^ 



./2 



-^0 



-0/2 

12/(^^2 + an) 



(55) 
(56) 

(57) 



If / > 0, then the matrix in (53 ) satisfies >- 0, so that the spectrum of J is purely 
imaginary. If, in addition, > 0, then det(JA^ — 0/2/2) 7^ and the first of the equations 
([55]) implies that = (in particular, k = 0) for the steady-state mean values of the 



system observables. Moreover, in this case, the matrix A in (57 1 is Hurwitz, thus leading 
to an admissible solution S := {crjk)i^j,k^2 (with S + iJ/2 :^ 0) of the second equation 
in ( [55] ) to be found from 

AS + SA^ + C - Ufau 



ail 
ail 2(712 



(58) 



where 



A :-- 



-0/2 
-^0 



(59) 



is a constant Hurwitz matrix. Due to the term, which is quadratic in S, the structure of 
(58 ) resembles that of an algebraic Riccati equation. In view of ( [59] ), the matrix algebraic 
equation (58), whose left-hand side is a real symmetric (2 x 2)-matrix, is equivalent to 
three equations 



2(712 - 00-11 + Cll 
-WgO-ii - 0(7i2 + (722 + C12 - 12/o-n 
-2Wo(Tl2 - 0(722 + C22 - 24/(Tii7i2 



0, 

0, 




(60) 
(61) 
(62) 



with three unknowns an, ai2, 022- By expressing cri2 and 0-22 from (60) and (62) in terms 

of o"ii as 



0-12 
0-22 



(0(711 - Cll)/2, 

(C22 - 2(u;2 + i2/aii)(7i2)/0 

((^0 + 12/(Tii)(cii - 0(7ii) + C22)/0, 



(63) 



(64) 



and substituting the representations into ( [61] ), it follows that an satisfies a quadratic equa- 
tion: 



24/a2, + {2ujI + 0V2 - 12/cii/0)(7ii - {{ujI + 0V2)cii + C22 + 0ci2)/0 = 0. (65) 
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Since the matrix (56) satisfies C ^ 0, so that C11C22 > c 
geometric mean inequality, 0^cii/2 + C22 ^ 0a/2ciiC22 > 



12' 



then, by the arithmetic - 



free term of the quadratic polynomial in (|65|) satisfies 

{V2 



a/20|ci2|. 

0V2)cii 



and hence, the 

f C22 + 0Ci2 > 



l)0|ci2| ^ 0. Therefore, in view of the assumption that / > 0, the polynomial 
has two real roots with opposite signs, of which the positive root an makes the matrix 



A in (57 1 Hurwitz and yields a unique admissible solution S for (58 1 whose other entries 
are computed through ( [63] ) and (|64]). As a numerical example, suppose the open quantum 
Duffing oscillator is driven by a four-dimensional quantum Wiener process, and 



Wo 



0.9026, B 



0.4853 
-0.5955 



-0.1497 
-0.4348 



-0.0793 
1.5352 



-0.6065 
-1.3474 



(66) 



Here, (54) is satisfied with (p = 0.6357. The behavior of the entries of the matrix S, com- 



puted by using the self-consistent quantum Gaussian linearization procedure for a range 
of nonnegative values of the anharmonicity parameter /, are shown in Fig. [T] The results 




Figure 1: The entries ajk of S, computed as an approximation to the real part of the 



steady-state quantum covariance matrix of the open quantum Duffing oscillator with (66) 
via the self-consistent Gaussian linearization for the range ^ / ^ 1 of the anharmonic- 
ity strength parameter in (52). The exact values of ajk in the purely harmonic case / = 
are marked by "o"s. 

of this approximation predict the decrease in the variance an of the position operator 
and the increase in the variance (T22 of the momentum operator as / increases. This is 



in qualitative agreement with the fact that the quartic term fq"^ in the Hamiltonian (52) 
with large positive / significantly "steepens" the walls of the potential well. In the case 
/ < 0, when the minimum of the potential V{q) at q = is only local, the properties of 
the linearization are more complicated and will be discussed elsewhere. 
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7 Conclusion 



We have proposed a quantum Gaussian linearization technique for a class of nonlinear 
open quantum systems with canonically commuting state variables governed by QSDEs 
with non-quadratic system Hamiltonians and bilinear coupling with the external heat bath. 
The approach is based on approximating the actual Hamiltonian by a quadratic function 
of the system observables in a mean square optimal fashion over a Gaussian reference 
quantum state. The optimal quadratic approximation of the Hamiltonian involves the 
inverse of a grade two special self-adjoint operator on real symmetric matrices, and we 
have described a method for its computation, more economical than that via vectorization 
of matrices. 

The resulting differential equations for the approximations of the mean and quan- 
tum covariance matrix of the system observables form a self-consistent set of equations 
which, despite their nonlinearity, produce a legitimate quantum covariance matrix. More- 
over, they preserve the CCRs of the system observables, thus making the effective linear 
dynamics physically realizable. We have demonstrated the approach for the quantum 
Duffing oscillator whose Hamiltonian is a quadro-quartic polynomial of the momentum 
and position operators. 

A time invariant version of the proposed technique involves nonlinear vector-matrix 
algebraic equations for which the existence/uniqeness of admissible solutions is, in gen- 
eral, a nontrivial problem to be tackled elsewhere. The error analysis and detailed dis- 
cussion of other properties of the quantum Gaussian linearization, and its applications 
to suboptimal filtering and control in nonUnear quantum systems, are also intended for 
further publications. 
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Appendices 

A Symmetries in the product moment of three observ- 
ables 

The following subsidiary lemma is used for establishing the symmetry of some of the 



mixed moments in ( 13 1, ( 14 1 and is provided here for completeness of exposition. 



Lemma 2 Suppose A, B, C are quantum observables on the underlying Hilbert space. 
Then the real part of their product moment satisfies 

ReE{ABC) = ReE{CBA). (Al) 
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Furthermore, if A and B (or B and C) commute canonically, then 

Re¥.{ABC) = ReE{BAC) (A2) 

(respectively, R.eE(ABC) = IieFj{ACB)). Finally, if each of the pairs {A,B) and 
{B, C) satisfies a CCR, then ReE{ABC) is invariant under arbitrary permutations of 
A, B, C. 



Proof. The first assertion ( Al I of the lemma follows from the identities 



E{ABC) = TiipABC) = Tt{CBAp) = TiipCBA) = E{CBA), (A3) 

where p is the density operator, and the self-adjointness of A, B, C is used. To prove 
the second statement of the lemma, we note that if the observables A and B commute 
canonically, that is, if \A, B] is a purely imaginary complex constant, then 

E{ABC) - E{BAC) = E{[A,B]C) = [A,B]EC 



is also purely imaginary by the realness of EC, thus implying (A2). The other case 
(when B and C commute canonically) is treated in a similar fashion, or alternatively, by 
reducing it to the previous case through repeatedly using the invariance of ReE{ABC) 
under swapping A with C. The third assertion of the lemma is established by combining 
the first two with the fact that the three transpositions A^C,A^B,B^C generate 
the group of all 6 possible permutations of A, B, C. ■ 



B Product moment of observables in a Gaussian quan- 
tum state 

Let C := (Ci)is£js£n be a vector of canonically commuting observables with a CCR matrix 
6 G A„, so that [C, C^] = iQ. We will compute the product moment E(^i x . . . x ^„) 
over a Gaussian quantum state in a different fashion from the traditional formulation 
and proof of Wick's theorem in terms of the annihilation and creation operators (|4]) and 
their normal ordering. To this end, by repeatedly using the Baker-Hausdorff formula 
gA+s ^ ^A^B-[A,B]/2 for operators A and B satisfying [A, [A, B]] = [B, [A, B]] = (see, 
for example, [8, pp. 128-129]), and the bilinearity of the commutator, it follows that 

n n 
k=l k=l 

for any u := {uj)i^j^n £ I^"- Here, ^5"fc=i A^ := Ai x . . . x An denotes the product 
of operators Ai, . . . , An, ordered "rightwards", with the order of multiplication being 
important in the noncommutative case. If the system is in a Gaussian quantum state, in 
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which ( has zero mean = and the quantum covariance matrix S := E(^^^), with 
ImS = 6/2, then the quantum characteristic function of ( is 



(B2) 



for any u G M", where S := Re^; cf. (|32]), (|33]). By combining (|B2]) with (|BT]), it follows 
that 



k=l 



Here, 



5* := {sjk)i<:jMn := S + iQ/2 



(B3) 



(B4) 



is a complex symmetric matrix, where := {9jk)i!ij,ksin is a real symmetric matrix which 
is defined in terms of the CCR matrix 6 as 



e 



kj 



e 



jk 



e 



jk, 



1 ^ j ^ k ^ n, 



(B5) 



and use is made of the property that all the diagonal entries of are zero. The product 
moment of the observables (^i, . . . , ^„ in the Gaussian quantum state is obtained from the 



function A in (B3 1 as 



(B6) 



fc=i 



Although the matrix S in (B4) is complex, its symmetry allows further analysis to be 
carried out similarly to the classical Gaussian case, as if S were the covariance matrix of 
a vector of jointly Gaussian real- valued random variables. Indeed, since A(— m) = A(m), 
all the odd order partial derivatives of A vanish at the origin, and hence, the product 
moment E^^^^^ in (B6 1 can only be nonzero if n is even, that is, if n = 2r for some 
positive integer r. In this case, a degree 2r homogeneous polynomial {—u^ Su/2Y /r\ in 
Ml, ... , U2r is the only term of the Taylor series expansion A(n) = Xlf^ol"^^'^^/^)^/^' 



of the right-hand side of ( |B3| ) which contributes to the partial derivative in ( |B6| ). More 
precisely, 



(9^1 .. . 9„2,.A(u)|^^g = coeff^ 



'-u^Su/2Y/r\ 



r!2' 



(B7) 



m=l 



where coeffi 
mial ui X . . . 



i...n2r(') denotes the coefficient of the polynomial associated with the mono- 
X U2r, and the sum is taken over a class of all possible (2r)! permutations 



, >, kr) of the integers 1, . . . , 2r. For any (ji, ki 



, jV , k^ ) 



G Tr, the prod- 
uct on the right-hand side of (B7 1 is invariant under r! permutations of the pairs (ji, ki), 
. . . , (jr, K) and with respect to r transpositions ji 



ki. 



Jr 



kr, with the latter 



invariance following from the symmetry of the complex matrix S defined by (B4)-(|B5[). 



The action of a group, generated by the r! pair permutations and r transpositions within 
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each of the pairs, splits into (2r)!/(r!2'') = (2r — 1)!! equivalence classes. Each of 
these classes is the orbit of the group passing through one of the regular permutations 
(ji, ki, . . . kr) of the integers 1, . . . , 2r which satisfies ji < ki, . . . , jr < kj. and 
ji < 32 < ■ ■ ■ < jr-i < jr (see Section |4]). Therefore, the sum in ( |B7[ ) reduces to that 
over the class Vr of regular permutations (as representatives of the equivalence classes) 



as 



dm ■ ■ ■ du,Mu)\u-. 



xn 



^ jm km ' 



(B8) 



Now, since Sjk - 
for any (ji,fci, 
follows that 



Sjk for all j ^ kin view of (B4)-(B5 1, then 11™=! ^ 



,jr, K) £ 'Pr- Hence, by combining (B6) (with n = 2r) and (B8), it 



2r 



k=l 



(jl,fcl,...,>,fcr)e'Pr ™=1 



(B9) 



which is what constitutes Wick's theorem. In the case = of pairwise commuting 
observables, when all the cross-covariances are symmetric, that is, Sjk = Skj, the relation 
( |B9[ ) reproduces Isserlis' theorem (W\ for jointly Gaussian random variables, which, due 
to the symmetry, is usually formulated in terms of pair partitions {{ji, ki}, . . . , {jr, K}} 
of the integers 1, . . . , 2r. However, in the noncommutative quantum case, when 6 7^ 
and Skj = 'sjk, the conditions ji < ki, jr < K in the definition of regular permuta- 



tions become essential for (B9 1 



C Proof of Lemma [11 



In the Gaussian quantum state, the expectation on the left-hand side of (20) can be com- 
puted for any matrix i? G §„ as 



j,k,£,m=l 
n 

= r jkV lm{SjkSim + SjiSkm + SjmSkl) 

j,k,£,m=l 

= {TT{RS)y + 2Tt{RS^RS) = {R, E)^ + 2Tr(i?(Ei?E + ei?0/4)) 

= {r,j:{j:,r) + 2K{R)), (ci) 

where we have used the symmetry of R and S and the antisymmetry of 9 (by which 
Tr(i?Ei?e) = 0), and also the representation (Issl) of the operator K{R) := Re{S'^RS). 



Comparison of the right-hand sides of (20) and (CI ) leads to (BTJ). We will now prove that 



S >~ entails K y 0. Since the property S y for the complex matrix (33 ) is equivalent 
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to the condition 



s 

0/2 



-0/2 
S 



y for the real matrices S and 0, it is also equivalent to 



that S :^ and that the real antisymmetric matrix 

E := S-i/2@5.-i/2/2 

is contractive in the sense that its operator norm ||S||oo 



(C2) 
(to be distin- 



guished from the Frobenius norm 
satisfies 



of matrices, with Amax(-) the largest eigenvalue) 



< 1. 



(C3) 



Since E is antisymmetric, its operator norm coincides with the spectral radius: ||S 
r(S). By bijectively transforming the matrix R into another real symmetric matrix 

R := V^rV^ 



(C4) 



(recall that S >- 0, which ensures the existence of a real positive definite symmetric matrix 
square root and combining (38) with (C2), it follows that 



{R, K{R)) = {R,R + ERE) = \\Rf + {R, ERE) ^ (1 



I ' — ^ II oo) 



li?ll 



(C5) 



Here, we have also used the Cauchy-Bunyakovsky-Schwarz inequality (applied to the 
Frobenius inner product of matrices) and the property that neither ||Hi?|| nor ||-RH|| ex- 
ceeds ||H||oo||-R||, by which 



{R,ERE) 



{ER,RE) ^ ^ 



|2 II D||2 
looll^ll • 



Indeed, \\ERf = Tt{RE'^ER) ^ A„ 



)Tr(i?2 
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I oo 



R\ 



and the inequahty 



^ ||S||oo ||-R|| is verified in a similar fashion. In view of (|C3j) and arbitrariness of 
the matrix i? G §„ in ( |C4[ ), the lower bound ( |C5| ) shows that the condition S y does 
entail K y 0. 



D Computing the inverse of the operator K in p8| ) 

Throughout this section, a short-hand notation I71, 5i \ . . . \ = Xlfc=il7A:; ^kl will 

be utilized for a special linear operator of grade r which acts on a matrix X as 

r 

|7i,(5i|...|7.,5.1(X):=^7fe^4, 

fc=i 

where 71, . . . , 7^, 5r are given appropriately dimensioned real matrices If for every 
k = 1, . . . ,r, the matrices 7/,., 5k are either both symmetric or both antisymmetric, I71, 5i \ 
. . . I 7r,5rl is a self-adjoint operator whose properties are studied in [30, Section 7]. 
Moreover, if for every k, the matrix 7^ has the same order and is either symmetric or 



'Such operator structure resembles the Kraus form of quantum operations 11211 pp. 360-373]. 
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antisymmetric, then |7i,7i | • • • | 7r,7rl is a self-adjoint operator on an appropriate 
subspace of real symmetric matrices. By using the identity I7, (5| |cr, r| = f'ja, t6J 
for the composition (interchangeably, product) of special operators of grade one, and 



recalling ( |C2| ), the operator K from ( [38[ ) can be factorized as 

K=lV^,V^miJ I S,Sl|^,v^l. (Dl) 

Here, = a/|S, S| is a positive definite special self-adjoint operator of grade 

one on the space S„ which was employed in \ca\ to carry out the transformation R h-^ R. 
The latter operator is straightforwardly invertible, with |v^, v^l"-^ = 
in view of S :^ 0. Therefore, ( |D1[ ) reduces the computation of the inverse operator 



(D2) 



to that of |/, / I S, Since the matrix S, defined by (C2), is real antisymmetric, 

it has purely imaginary spectrum ±icji, . . . , ±iuj^ (with all coi, ... ,uj^ real and, without 
loss of generality, nonnegative) and is orthogonally block-diagonalizable in the sense that 

E = U{3® 13)U^. (D3) 

Here, U E M"^" is an orthogonal matrix whose columns are formed from the real and 
imaginary parts of u eigenvectors of S associated with the eigenvalues iui, . . . , ico^. Also, 
the matrix J is given by Q, and 



15 := diag {uk] 



(D4) 



is a diagonal matrix with ui, . 
product of matrices. 



, over the main diagonal, so that, with (g) the Kronecker 
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U 
-U 



(D5) 



is a two-diagonal real antisymmetric matrix. Here, we prefer to utilize real matrices in 
order to avoid extensions of linear operators to spaces of complex matrices. The orthogo- 
nality of the matrix U implies that f/^| is a unitary operator on the space §„. Indeed, 
since its adjoint is fU, = fU^, Uj, then 



lu, u^^iu, uq = lU'U, U'U} = II, 11= I 



is the identity operator on §„. From (D3 ) and the unitarity of |f/, f/^|, it follows that 

|f/^f/l|/,/|H,Sir,f/l=X+Z, 

where 



(D6) 



Z := |J ® O, J ® 01 (D7) 
is a grade one self-adjoint operator on §„, whose operator norm coincides with its spectral 



radius and is computed in terms of ( |D4[ )-( p5) ) as 

r{Z) = r(J®0)2 = 



max Uu. 
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Therefore, since the matrix S is contractive in view of ( |C3| ), so that ^ Wi, . . . , < 1, 
then so also is the operator thus ensuring the invertibility of the operator X + Z. Hence, 
in view of the unitarity of U\ in (D6), the inverse operator K^^ in (D2) can be 
computed as 

K-^ = |s-v^s-l/l|f/,f/^(x+z)-l|f/^f/lls-v^s-l/l 

= |f/,f/l(X+Z)-i|f/^f/l. (D8) 

U := T.-^/'^U (D9) 

is a nonsingular matrix which satisfies Tj'^QU /2 = [/( J ® O) and is, therefor e, re lated to 
the eigenvectors of the matrix S^^0/2 = E~^/^Sv^, isospectral to S from (C2). Now, 
to compute the inverse operator on the right-hand side of (|D8[), we factorize it as 



Here, 



where 



and 



D := diag (4) 

IsSfcsSn 





u"^ 



(DIO) 
(DID 
(D12) 



is a diagonal matrix with diagonal entries = d^j^y = tul for k = 1, . . . , z/. Here, we 
have used the property = —I2 for the matrix ^ and the power identities I7, ^l''" = 
17'^, 5''J and (7 5)'' = 7^ (J'^ for square matrices 7, 5 and nonnegative integers k. It 
follows from (Dll )-(D12) that the inverse operator on the right-hand side of (DIO) can 
be expanded into an absolutely convergent operator power series as 



(D13) 



Due to the diagonal structure of the matrix D from (D12), the image (X — Z'^)^^{Y) of a 
matrix Y := {yjk)i^j,k^n G §n under the operator (D13) is a real symmetric matrix with 
entries 

{{I - Z')-\Y)),k = Y.^]y,kdi = y,k/il - d,du) 



for all 1 ^ j,k ^ n. Hence, upon splitting the matrix Y :- 
{u X z/)-blocks Yjk, its image takes the form 



Y21 Y22 



{X-Z')-\Y) 



LQYii LQ Y12 
LQY21 LQ Y22 



(I2 0L)QY. 



into four 



(D14) 



Here, denotes the Hadamard product of matrices, the matrix L EE>yis defined by 

L:={l/{l-ujy,))i^,,k^y, (D15) 
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and Ir denotes the (r x r)-matrix of ones. To represent the other factor I — Z on the 



right-hand side of (DIO), we note that (|D4[) and (D5) allow the image of the matrix Y 



under the operator Z in ( |D7[ ) to be computed as 

Z{Y) = 



u 
-U 



Y 



U 
-U 



-M Y22 
MQYu 



MQY21 
-M Yii 



-UY22U UY21U 
UY12U -UYnU 



[I2 ® M) ((J ® QY{3 ® /,)),(D16) 



where 



M :-- 



[(jJjUJk)l^j,k^u — , 



(D17) 



Finally, by assembling ( |D8[ ), ( |D10[ ), ( |D14[ ), ( |D16[ ) together, it follows that the inverse 
operator K^^ can be computed for the matrix F G §„ in (39) as 



K~\T) = f/((l2 ® L) (F - (I2 M) ((J /,)F(J Q)))U^ 



(D18) 



with F := U^TU. Here, the matrices U, L, M are related by (|D9|), ( |D15| ), ( |D17| ) with 
the eigenvectors and eigenvalues of the matrix S from (|C2|). It now remains to note that 



the positiveness (40) of K^^ follows from (D18) which represents this operator as the 
composition of positive operators |[7'^,[7|, |J'^0/^, J0/^|,X+ (l20M)0, (l20L)0 
and |f/, f/^|. Here, we have used the positive semi-definiteness of the matrices 1^, M 
and 



L = l^ + y^ M0 0M 



k^l 2k times 

and the fact that the Kronecker and Hadamard products of positive semi-definite matrices 
are also positive semi-definite; see the Schur product theorem [i9i Theorem 7.5.3 on p. 
458]. 



E Quadratic approximation of the Hamiltonian for the 
quantum Duffing oscillator 



The mean square optimal quadratic approximation of the Hamiltonian iJ in (52) reduces 
to the approximation of the quartic term by a quadratic function of q and p. The latter 
problem reduces to the quadratic approximation of the non-quadratic part of 



A:= Akx^ + x\ 
^ V ' 

non— quadratic 



(El) 



in the reference quantum state, with x and w the centered position and momentum oper- 
ators: 

w 



X 



q- K, 



K, := Eg, 



zu := p — Ep, 



(E2) 
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where the centered vector ^ of system observables is defined in conformance with ([T2]). 



We will now compute the contribution of A from (El) to the quadratic approximation 



of the Hamiltonian ( p2| ) through the vector e and the matrix F from ( [13] ) in a Gaussian 

ReEirJ^e) = ReE(x'a^) - Sa^.S. (E3) 



reference quantum state: 



r 



Here, in accordance with ( [T2| ) and ( |E2[ ), and by applying Wick's theorem (see Appendix|B]) 

to E(x^) = 3(E(x2))2, 

?7 := A - EA = A - 3ai\ (E4) 

is obtained by centering the non-quadratic part A, with an := E(x^) the variance of 
the position operator, and use is made of the relations E(x'^^) = and E(x^^^^) = 
which follow from the property that the odd order central moments vanish in a Gaussian 
quantum state. Throughout this section, the "hat" symbol marks the quantities associated 
with the quadratic approximation of A, such as 'e, T, fj in (E3) and (E4). They will be 
multiplied by / and combined with the remaining quadratic part of the Hamiltonian in 
([52]). Application of Wick's theorem to the mixed moments in ([E3]) yields 



3E(x')E(xO = 3a 



11 



3(Tii 



3(E(x'))'E(e,efc) + 12E(x')E(xe,)E(xe 
3o"iiSjfc + 12criiSijSifc 



0-11 
Sl2 



for all 1 ^ j, A; ^ 2, where Sjk denote the entries of the quantum covariance matrix 



S 



Sll Si2 




' E(x2) 






0-11 0-12+^2 


S21 S22 






E(tI72) 




^■12 — V2 <722 



(E5) 
(E6) 

(E7) 



obtained from (]3]), <\33^, with S : 

is Eix'^e: 



(o"jfe)is;j,fc^2 G A vector-matrix 



'orm of (E6) 



Sl2 



condition S Q of Theorem [2] if and only if an > 0, (T22 > and 

det S = criiO-22 - o-?2 > 1/4, 



[ Sll S12 ] . The matrix S from (E7) satisfies the 

(E8) 



which is stronger than S :^ 0. By combining ([E5[)-([E7[), it follows that ( [E3[ ) take the form 

^= 12k(Tii(J, 



r = Sail Aaa 




1 



(E9) 



where a 



o-ll 
0-12 



denotes the first column of the symmetric matrix S. Since E a is 



the first column of I2, substitution of ( [E9[ ) into ( ]39] ) yields 

/So = 12fi;criiS^V = 12Kaii 



(ElO) 
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where k is the mean value of the position operator from (E2). We will now compute the 
matrix in (39) by employing (D18) of Appendix [D| Since the matrix J from (|3]) spans 
the space A2, the eigenvalues iiwi of the real antisymmetric (2 x 2)-matrix 



(Ell) 



from (C2) and (D3 1 (with U G M^^^ an orthogonal matrix]^ being specified by the eigen- 



vectors) are given by 



1 



2 Vdet S 



< 1, 



(E12) 



where the inequality follows (E8). In the example being considered, the matrices L and 
M from (D15) and (D17) become positive scalars: L = 1/(1 — ujf) and M = cof, so 
that the action of each of the operators (I2 ® L)Q and (I2 ® M)0 on a (2 x 2)-matrix 



is equivalent to an appropriate scaling of the matrix. This allows (D18 1, in application to 
([39]), to be simplified as 



Ro 



1 



jrj 



(2detS)2 



(E13) 



with 



-1 



Sail 4 



1 




^-1 




1 



^-1 



(E14) 



in view of (|E9|). Here, T := U'^TU and U := T.'^/'^U in acco rdanc e with (|D9|), so that 
UU'^ = in view of the orthogonality of the matrix U from (El 1 1. Also, use has been 
made of the identity Y^^ = —3Y3/ detV which holds for any nonsingular symmetric 



matrix Y of order two. Substitution of (E9), (E12), (E14) into ( E13| ) yields 



Ro = I2au 



1 




(E15) 



In view of (ElO) and (El 5 1, the mean square optimal quadratic approximation (3^ + 
Roil 2 of the operator A in (El I (with the additive constant terms being omitted) does 
not depend on the momentum operator. Now, by multiplying the quadratic approximation 
of A by / and combining the result with the remaining quadratic part of the Hamilto- 



nian (52), it follows that its mean square optimal quadratic approximation in the Gaussian 
quantum state is 



H 



^-^^ + fiAn'x + Q^\' + + eRM2) + (*) 



^The matrix U is also symplectic with the structure matrix J in the sense that C/J[/^ = J. 
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where (*) assembles the additive constant terms which are irrelevant for the dynamics of 
the quantum Duffing oscillator, and 



Ep 

Lo^ + 12/k2 
1 



+ fPo = 



a;2 + 12/(^2 ^ an) 
1 



(E16) 
(E17) 
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